suppressPackageStartupMessages({
library(SummarizedExperiment)
library(SEtools)
library(edgeR)
library(DT)
library(pheatmap)
library(plotly)
library(dplyr)
library(sva)
library(reshape2)
library(overlapper)
library(ggrepel)
})
source("Functions/EDC_Functions.R")
source("Functions/CriFormatted_OLD.R")
source("Functions/CriFormatted.R")
load("Data/TotalSumExp.RData")
load("Data/geneLengths.RData")
load("Data/AllSEcorrected.RData",verbose=T)## Loading objects:
## SEs
## DEAs
load("Data/brainspan.cortex.byAge.RData", verbose = T)## Loading objects:
## byAge
load("Data/iPSC.RData", verbose = T)## Loading objects:
## e
Organoids <- assays(SEs$chronic.org[,which(SEs$chronic.org$EXPO=="CNT")])$counts
Fetal2D <- assays(SEs$chronic.fetal[,which(SEs$chronic.fetal$EXPO=="CNT")])$counts
iPSC <- e[,grep("3391B",colnames(e))]
Fetal <- byAge
i <- intersect(rownames(Organoids),intersect(rownames(iPSC),rownames(Fetal2D)))
tot <- cbind(Organoids[i,],iPSC[i,],Fetal2D[i,])
tot <- filterGenes(tot)
en <- donorm(tot)
#put the rows in the same order
i <- intersect(row.names(Fetal),row.names(en))
en <- en[i,]
Fetal <- Fetal[i,]
#create the matrix of pearson correlation coefficients and plot
m <- matrix(0,nrow=ncol(Fetal),ncol=ncol(en))
for(i in 1:ncol(en)){ m[,i] <- apply(Fetal,2,FUN=function(x){ cor(log(en[,i]+1),log(x+1)) }) }
row.names(m) <- colnames(Fetal)
colnames(m) <- colnames(en)
m2 <- m[grep("pcw",row.names(m)),]
byheatmap(m2,cluster_row=F)org.chronic <- row.names(DEAs$chronic.org)[which((abs(DEAs$chronic.org$logFC.EXPO1X)>0.5 | abs(DEAs$chronic.org$logFC.EXPO1000X)>0.5) & DEAs$chronic.org$FDR<=0.05 & DEAs$chronic.org$logCPM > 0)]
MixN_org_chronic <- SEs$chronic.org[,which(SEs$chronic.org$EXPO=="CNT"|SEs$chronic.org$EXPO=="1X" | SEs$chronic.org$EXPO=="1000X")]
sehm(MixN_org_chronic, assayName = "corrected", org.chronic, do.scale = T, show_rownames = F, anno_columns = c("Line","EXPO2"), main="Chronic org DEGs")FDRTh <- 0.05
LogFCTh <- 0.5
MixN1X <- DEAs$chronic.org
MixN1X$genes <- row.names(MixN1X)
MixN1X$logFC <- (MixN1X$logFC.EXPO1X)
MixN1000X <- DEAs$chronic.org
MixN1000X$genes <- row.names(MixN1000X)
MixN1000X$logFC <- (MixN1X$logFC.EXPO1000X)
#Comp <- compareResultsFC(MixN1X, MixN1000X, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, logCPMth = 0, FCceil=8, title='Comparing MixN 1X to 1000X', geneLabel=TRUE)
#Comp$ScatterComp <- compareResultsFCNew(MixN1X, MixN1000X, FDRth=0.05, FCth=2^0.5, FDRceil=1e-10, logCPMth=0, FCceil=2.5,
title='Comparing MixN 1X to 1000X', geneLabel=TRUE, topLab=-16)
## Warning in cor.test.default(AllRes$logFC_A, AllRes$logFC_B, method = corMethod):
## Cannot compute exact p-value with ties
Comp$Scatter
ggsave(Comp$Scatter, filename='Data/ChronicOrganoids/Scatter.pdf', width=10, height=10)
ggsave(Comp$Scatter, filename='Data/ChronicOrganoids/Scatter.png', width=10, height=10)load("Data/ASD.RData", verbose = T)
## Loading objects:
## ASD
## SFARI
## SFARIgenes
## AutismKB
## DecipherNeuro
## MSSNG
## ASDsevere
## ID
## NeuroOmim
## iPSYCH
## GDI
## NeuropsychiatricDiseases
org.chronic3 <- row.names(DEAs$chronic.org)[which(DEAs$chronic.org$FDR<=0.05)]
controlNeg <- row.names(DEAs$chronic.org)[-which(DEAs$chronic.org$FDR<=0.05)]
random1 <- sample(row.names(DEAs$chronic.org), 10)
random2 <- sample(row.names(DEAs$chronic.org), 100)
random3 <- sample(row.names(DEAs$chronic.org), 1000)
chronic.orgUp <- row.names(DEAs$chronic.org)[which(((DEAs$chronic.org$logFC.EXPO1X)>0 & (DEAs$chronic.org$logFC.EXPO1000X)>0) & DEAs$chronic.org$FDR<=0.05)]
chronic.orgDown <- row.names(DEAs$chronic.org)[which(((DEAs$chronic.org$logFC.EXPO1X)<0 & (DEAs$chronic.org$logFC.EXPO1000X)<0) & DEAs$chronic.org$FDR<=0.05)]
MixN <- list(organoidsDEGs=org.chronic3, organoidsDEGsUp=chronic.orgUp, organoidsDEGsDown=chronic.orgDown,NonAffectedGenes=controlNeg, random1=random1,random2=random2,random3=random3 )
m <- overlapper::multintersect(ll = MixN, ll2 = ASD, universe = row.names(DEAs$chronic.fetal))
#overlapper::dotplot.multintersect(m,sizeRange = c(0,15))
org.chronicLong<-org.chronic3
save(org.chronic,org.chronicLong,file = "Data/DEGsOrganoidsChronic.RData")dotplot.multintersect.mod(m, sizeRange = c(0,15), th=0.05)
ggsave(dotplot.multintersect.mod(m, th=0.05), filename='Data/ChronicOrganoids/DotplotFetalNDD.pdf', width=10, height=6)m$prob <- round(m$prob, digits=3)
#DT::datatable(m$prob)m <- overlapper::multintersect(ll = MixN, ll2 = SFARIgenes, universe = row.names(DEAs$chronic.fetal))dotplot.multintersect.mod(m, th=0.05)
ggsave(dotplot.multintersect.mod(m, th=0.05), filename='Data/ChronicOrganoids/DotplotFetalSFARI.pdf', width=10, height=6)m <- overlapper::multintersect(ll = MixN, ll2 = NeuropsychiatricDiseases, universe = row.names(DEAs$chronic.fetal))dotplot.multintersect.mod(m, th=0.05)
ggsave(dotplot.multintersect.mod(m, th=0.05), filename='Data/ChronicOrganoids/DotplotFetalNPD.pdf', width=10, height=6)https://bioconductor.org/packages/release/bioc/html/GeneOverlap.html
library(GeneOverlap)
gom.obj <- newGOM(MixN,ASD, genome.size=length(row.names(DEAs$chronic.org)))
print(gom.obj)
## A GeneOverlapMatrix object:
## ###### Intersection ######
## SFARI AutismKB Decipher MSSNG iPSYCH NeuroOMIM ASDsevere ID
## organoidsDEGs 172 39 21 15 29 15 40 43
## organoidsDEGsUp 53 7 9 4 10 9 14 12
## organoidsDEGsDown 68 20 7 10 13 2 19 21
## NonAffectedGenes 673 178 134 39 73 65 167 144
## random1 0 0 0 0 0 0 0 0
## random2 4 0 2 0 0 1 0 0
## random3 56 11 12 3 8 4 13 18
## SingleCellASD NeuronASD AstrocyteASD NeuropsychiatricDiseases
## organoidsDEGs 109 26 8 226
## organoidsDEGsUp 24 5 0 76
## organoidsDEGsDown 57 17 7 87
## NonAffectedGenes 344 69 41 935
## random1 0 0 0 1
## random2 2 0 1 6
## random3 33 7 3 65
## ###### P-value ######
## SFARI AutismKB Decipher MSSNG iPSYCH
## organoidsDEGs 0.02231899 0.2882567 0.8344645 0.018552716 0.0006840381
## organoidsDEGsUp 0.14824495 0.9237283 0.3571956 0.268247645 0.0239544632
## organoidsDEGsDown 0.47977370 0.1786488 0.9380662 0.004489613 0.0260206988
## NonAffectedGenes 1.00000000 0.9958340 0.6417477 0.996755940 0.9997043028
## random1 1.00000000 1.0000000 1.0000000 1.000000000 1.0000000000
## random2 0.84602444 1.0000000 0.2570886 1.000000000 1.0000000000
## random3 0.66744915 0.8328923 0.2765499 0.663609205 0.2875722940
## NeuroOMIM ASDsevere ID SingleCellASD NeuronASD
## organoidsDEGs 1.0000000 0.4769800 0.00783293 0.0003390561 0.015528828
## organoidsDEGsUp 0.9554060 0.3202309 0.20909608 0.5742959087 0.601455098
## organoidsDEGsDown 1.0000000 0.4070359 0.02858962 0.0003726425 0.001505528
## NonAffectedGenes 1.0000000 1.0000000 0.99976073 1.0000000000 0.999999983
## random1 1.0000000 1.0000000 1.00000000 1.0000000000 1.000000000
## random2 0.8465324 1.0000000 1.00000000 0.8269194222 1.000000000
## random3 0.9999934 0.7779795 0.04787005 0.4122299829 0.505493712
## AstrocyteASD NeuropsychiatricDiseases
## organoidsDEGs 0.7377828 0.0007054163
## organoidsDEGsUp 1.0000000 0.0061456247
## organoidsDEGsDown 0.1258210 0.3705189362
## NonAffectedGenes 0.9994944 0.9999999973
## random1 1.0000000 0.5324600843
## random2 0.3089789 0.7491843066
## random3 0.7203668 0.8626548875
## ###### Odds Ratio ######
## SFARI AutismKB Decipher MSSNG iPSYCH NeuroOMIM
## organoidsDEGs 1.1983919 1.1192904 0.8166280 2.0367785 2.1644293 0.27878640
## organoidsDEGsUp 1.1784266 0.6274643 1.1891777 1.5648136 2.1796504 0.60720421
## organoidsDEGsDown 1.0129806 1.2786952 0.6036833 2.9608197 1.9457765 0.08673027
## NonAffectedGenes 0.4071181 0.6529367 0.9512470 0.4485827 0.4620156 0.04618727
## random1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000000
## random2 0.6664929 0.0000000 2.0737493 0.0000000 0.0000000 0.53402587
## random3 0.9478374 0.7737971 1.2445912 0.8831531 1.3057397 0.20237707
## ASDsevere ID SingleCellASD NeuronASD AstrocyteASD
## organoidsDEGs 1.0220948 1.5724179 1.4816506 1.7034684 0.8325990
## organoidsDEGsUp 1.1746139 1.3327842 0.9767273 0.9571983 0.0000000
## organoidsDEGsDown 1.0816286 1.6385936 1.6902996 2.4694903 1.7542714
## NonAffectedGenes 0.3550056 0.5493021 0.3599084 0.3151439 0.3967105
## random1 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## random2 0.0000000 0.0000000 0.6271536 0.0000000 2.7654528
## random3 0.8307804 1.5948446 1.0546567 1.0513365 0.8054194
## NeuropsychiatricDiseases
## organoidsDEGs 1.2900641
## organoidsDEGsUp 1.3919764
## organoidsDEGsDown 1.0445669
## NonAffectedGenes 0.6423804
## random1 1.4073308
## random2 0.8073267
## random3 0.8734341
## ###### Jaccard Index ######
## SFARI AutismKB Decipher MSSNG iPSYCH
## organoidsDEGs 0.051667167 0.014275256 0.007829978 0.005807201 0.011085627
## organoidsDEGsUp 0.031454006 0.006993007 0.009667025 0.004813478 0.011467890
## organoidsDEGsDown 0.033415233 0.014781966 0.005392912 0.008403361 0.010534846
## NonAffectedGenes 0.047879909 0.012878952 0.009712256 0.002828752 0.005289855
## random1 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000
## random2 0.003795066 0.000000000 0.007751938 0.000000000 0.000000000
## random3 0.029442692 0.009038620 0.010452962 0.002851711 0.007312614
## NeuroOMIM ASDsevere ID SingleCellASD NeuronASD
## organoidsDEGs 0.005300353 0.01451906 0.01597325 0.036986766 0.009900990
## organoidsDEGsUp 0.008387698 0.01375246 0.01250000 0.018912530 0.005656109
## organoidsDEGsDown 0.001384083 0.01378810 0.01595745 0.035602748 0.013742926
## NonAffectedGenes 0.004640206 0.01205254 0.01042044 0.024677188 0.004996018
## random1 0.000000000 0.00000000 0.00000000 0.000000000 0.000000000
## random2 0.002493766 0.00000000 0.00000000 0.003273322 0.000000000
## random3 0.003081664 0.01049233 0.01533220 0.022297297 0.006352087
## AstrocyteASD NeuropsychiatricDiseases
## organoidsDEGs 0.003082852 0.0643691256
## organoidsDEGsUp 0.000000000 0.0400421496
## organoidsDEGsDown 0.005843072 0.0386323268
## NonAffectedGenes 0.002973169 0.0666429081
## random1 0.000000000 0.0008312552
## random2 0.006289308 0.0046583851
## random3 0.002838221 0.0305307656
drawHeatmap(gom.obj)gom.obj <- newGOM(MixN,SFARIgenes, genome.size=length(row.names(DEAs$chronic.org)))
drawHeatmap(gom.obj)gom.obj <- newGOM(MixN,NeuropsychiatricDiseases, genome.size=length(row.names(DEAs$chronic.org)))
drawHeatmap(gom.obj)geneStripPairEDCMix(SE = MixN_org_chronic,GeneSet =intersect(org.chronic, union(union(SFARIgenes$score1,SFARIgenes$score2),SFARIgenes$score3)), printExp = FALSE, SampleColors = "Default")
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.ymin` is deprecated. Use `fun.min` instead.
## Warning: `fun.ymax` is deprecated. Use `fun.max` instead.PVals <- DEAs$chronic.org[intersect(org.chronic, unlist(SFARIgenes)),]
PVals <- PVals[complete.cases(PVals), ]
PVals
## logFC.EXPO1X logFC.EXPO1000X logCPM F PValue
## SCN9A 0.16332099 -0.73394503 3.6238927 28.679243 4.486322e-08
## NEFL 0.17224825 -0.67496037 7.2905518 26.587394 1.004741e-07
## CACNA1E -0.22395439 -0.68438361 4.9734221 21.775175 7.506001e-07
## SCN2A -0.23500980 -0.75882170 4.4730437 21.537869 8.340935e-07
## PCDH9 -0.35329648 -0.51191575 6.2629341 18.479058 3.451823e-06
## CDH9 -0.63189859 -1.50276593 0.4498228 17.344314 6.030875e-06
## SPARCL1 0.28517770 -0.59496470 5.7544922 15.403606 1.636156e-05
## TH 0.31734183 -0.56675217 2.3380022 14.953642 2.079647e-05
## NEGR1 -0.04352262 -0.54219055 4.4452160 14.906921 2.132497e-05
## SLC6A1 -0.21067534 -0.57822638 3.0876260 13.077681 5.868259e-05
## NRXN3 -0.50939333 0.21190259 5.1255620 12.755676 7.057116e-05
## KCNMA1 -0.65819131 -0.33457106 2.9964461 12.306480 9.158739e-05
## GRM8 0.59891403 0.30232488 5.4549395 11.619416 1.375170e-04
## CACNB2 -0.25652132 -0.51771213 3.2092017 11.613580 1.379983e-04
## DPYD -0.41965092 -0.59127083 3.0251296 11.342770 1.624122e-04
## LMX1B 0.22470186 -0.62563738 4.3620879 11.132867 1.844648e-04
## KDM6B 0.28440194 0.54084969 5.1366501 10.480879 2.756171e-04
## WNT1 0.42357343 -0.73059040 1.7103686 10.118513 3.459453e-04
## FABP3 -0.54794627 -0.50839479 3.3966961 9.978136 3.780888e-04
## CNTN5 0.05897304 -0.60625916 1.6670097 8.563565 9.503601e-04
## DRD2 0.17067356 -0.65784667 1.0232263 8.552201 9.576166e-04
## KIRREL3 -0.06077077 -0.51471524 2.7742917 8.078676 1.318723e-03
## PCDH11X -0.32805473 -0.64062334 4.6909851 7.525175 1.931552e-03
## GRM7 0.05663132 -0.61625257 2.1733713 7.382859 2.133643e-03
## SLC4A10 0.10963873 0.57241339 5.2768430 7.242564 2.354874e-03
## ACE 0.50518870 -0.00501186 1.4608377 7.221730 2.389741e-03
## FABP7 -0.58477173 -0.21714537 4.9423324 6.660712 3.567096e-03
## THBS1 -0.55858202 -0.12754071 3.7180289 6.643538 3.611637e-03
## FEZF2 0.33938486 0.89684384 2.2713893 6.558978 3.839691e-03
## ADAMTS18 0.30426384 0.60810219 3.1172417 6.255539 4.791842e-03
## SAMD11 0.61769439 0.16814633 2.5613572 5.718011 7.145064e-03
## CACNA1I -0.56456711 -0.82961573 1.3985198 5.689791 7.298381e-03
## FOXP2 -0.57048773 0.03921951 5.9517196 5.616600 7.712461e-03
## FDR
## SCN9A 0.0001156523
## NEFL 0.0001170810
## CACNA1E 0.0003182873
## SCN2A 0.0003239857
## PCDH9 0.0006120982
## CDH9 0.0008920245
## SPARCL1 0.0014428240
## TH 0.0016448800
## NEGR1 0.0016487940
## SLC6A1 0.0029098719
## NRXN3 0.0031645244
## KCNMA1 0.0036541645
## GRM8 0.0047692248
## CACNB2 0.0047697124
## DPYD 0.0052581297
## LMX1B 0.0057095415
## KDM6B 0.0070921402
## WNT1 0.0081205050
## FABP3 0.0084843741
## CNTN5 0.0146529435
## DRD2 0.0147116784
## KIRREL3 0.0175908790
## PCDH11X 0.0222464372
## GRM7 0.0234873475
## SLC4A10 0.0246898511
## ACE 0.0249591803
## FABP7 0.0313881338
## THBS1 0.0316627215
## FEZF2 0.0327906209
## ADAMTS18 0.0371020919
## SAMD11 0.0473454810
## CACNA1I 0.0477600401
## FOXP2 0.0495748998load("Data/DEGsFetalChronic.RData", verbose = T)
## Loading objects:
## fetal.chronic
## fetal.chronicLong
geneStripPairEDCMix(SE = MixN_org_chronic,GeneSet = intersect(fetal.chronic, unlist(SFARIgenes)), printExp = FALSE,SampleColors = "Default")
## [1] "Expression values are not available for the following genes: CD38"
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.ymin` is deprecated. Use `fun.min` instead.
## Warning: `fun.ymax` is deprecated. Use `fun.max` instead.PVals <- DEAs$chronic.org[intersect(fetal.chronic, unlist(SFARIgenes)),]
PVals <- PVals[complete.cases(PVals), ]
PVals
## logFC.EXPO1X logFC.EXPO1000X logCPM F PValue FDR
## RANBP17 -0.1626618 -0.2152368 5.277389 4.470447 0.018742025 0.08479478
## UPF3B 0.0328277 -0.2586830 3.815079 2.967897 0.064606056 0.17889598
## KCNJ2 0.4778245 0.1451159 6.134523 5.280253 0.009961453 0.05816433load("Data/DEGsFetalAcute.RData", verbose = T)
## Loading objects:
## fet.acute
geneStripPairEDCMix(SE = MixN_org_chronic,GeneSet = intersect(fet.acute,org.chronic), printExp = FALSE,SampleColors = "Default")
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.ymin` is deprecated. Use `fun.min` instead.
## Warning: `fun.ymax` is deprecated. Use `fun.max` instead.PVals <- DEAs$chronic.org[intersect(fet.acute,org.chronic),]
PVals <- PVals[complete.cases(PVals), ]
PVals
## logFC.EXPO1X logFC.EXPO1000X logCPM F PValue FDR
## NEFL 0.1722482 -0.6749604 7.290552 26.58739 1.004741e-07 0.000117081geneStripPairEDCMix(SE = MixN_org_chronic,GeneSet = c("CLSTN2", "EPHB2"), printExp = FALSE,SampleColors = "Default")
## Warning: `fun.y` is deprecated. Use `fun` instead.
## Warning: `fun.ymin` is deprecated. Use `fun.min` instead.
## Warning: `fun.ymax` is deprecated. Use `fun.max` instead.PVals <- DEAs$chronic.org[c("CLSTN2", "EPHB2"),]
PVals <- PVals[complete.cases(PVals), ]
PVals
## logFC.EXPO1X logFC.EXPO1000X logCPM F PValue FDR
## CLSTN2 0.214876852 0.09619874 5.239390 1.447981 0.24891393 0.4130782
## EPHB2 0.006988673 0.16640719 7.820478 3.802717 0.03213182 0.1164109load("Data/DEGsFetalChronic.RData", verbose = T)
## Loading objects:
## fetal.chronic
## fetal.chronicLong
load("Data/DEGsOrganoidsChronic.RData", verbose = T)
## Loading objects:
## org.chronic
## org.chronicLong
MixN <- list(OrganoidsDEGs=org.chronic, FetalAcuteDEGs=fet.acute, FetalChronicDEGs=fetal.chronic)
m<-multintersect(ll = MixN, universe = row.names(filterGenes(Total)))
dotplot.multintersect.mod(m, th=0.05)
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).m<-multintersect(ll = MixN, ll2 = ASD, universe = row.names(filterGenes(Total)))
dotplot.multintersect.mod(m, th=0.05)m<-multintersect(ll = MixN, ll2 = SFARIgenes, universe = row.names(filterGenes(Total)))
dotplot.multintersect.mod(m, th=0.05)m<-multintersect(ll = MixN, ll2 = NeuropsychiatricDiseases, universe = row.names(filterGenes(Total)))
dotplot.multintersect.mod(m, th=0.05)load("Data/DEGsFetalChronic.RData", verbose = T)
## Loading objects:
## fetal.chronic
## fetal.chronicLong
load("Data/DEGsOrganoidsChronic.RData", verbose = T)
## Loading objects:
## org.chronic
## org.chronicLong
MixN <- list(OrganoidsDEGs=org.chronicLong, FetalAcuteDEGs=fet.acute, FetalChronicDEGs=fetal.chronicLong)
m<-multintersect(ll = MixN, universe = row.names(filterGenes(Total)))
dotplot.multintersect.mod(m, th=0.05)
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).m<-multintersect(ll = MixN, ll2 = ASD, universe = row.names(filterGenes(Total)))
dotplot.multintersect.mod(m, th=0.05)m<-multintersect(ll = MixN, ll2 = SFARIgenes, universe = row.names(filterGenes(Total)))
dotplot.multintersect.mod(m, th=0.05)m<-multintersect(ll = MixN, ll2 = NeuropsychiatricDiseases, universe = row.names(filterGenes(Total)))
dotplot.multintersect.mod(m, th=0.05)